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If a deuterated molecule containing strong intramolecular hydrogen bonds is placed in a 
hydrogenated solvent it may preferentially exchange deuterium for hydrogen. This prefer¬ 
ence is due to the difference between the vibrational zero-point energy for hydrogen and deu¬ 
terium. It is found that the associated fractionation factor <1) is correlated with the strength of 
the intramolecular hydrogen bonds. This correlation has been used to determine the length 
of the H-bonds (donor-acceptor separation) in a diverse range of enzymes and has been ar¬ 
gued to support the existence of short low-barrier H-bonds. Starting with a potential energy 
surface based on a simple diabatic state model for H-bonds we calculate $ as a function of 
the proton donor-acceptor distance R. For numerical results, we use a parameterization of 
the model for symmetric 0-H- • • O bonds [Chem. Phys. Lett. 535, 196 (2012)]. We con¬ 
sider the relative contributions of the 0-H stretch vibration, 0-H bend vibrations (both in 
plane and out of plane), tunneling splitting effects at finite temperature, and the secondary 
geometric isotope effect. We compare our total as a function of R with NMR experi¬ 
mental results for enzymes, and in particular with an empirical parametrization ^{R), used 
previously to determine bond lengths. 
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I. INTRODUCTION 


The issue of low-barrier hydrogen bonds in proteins and whether they play any funetional role, 


partieularly in enzyme eatalysis, is eontroversial llll-llflh. Identifying sueh short hydrogen bonds, 


eharaeterised by a donor-aeeeptor distanee of i? ~ 2.45—2.65 A, is not eompletely straightforward 
111 in . In protein X-ray erystallography, the standard errors in inter-atomie distanees are about 10 
to 30 per eent of the resolution. Henee, for an X-ray strueture with 2.0 Aresolution, the standard 
errors in the distanees are ±(0.2 —0.6) A. This uneertainty makes it diffieult to distinguish between 
short strong bonds and the more eommon weak long bonds, with R > 2.8 A [I12I] . NMR provides 
an alternative method of bond length determination via the ehemieal shift. An independent 
NMR “ruler” involves isotopie fraetionation, where one measures how mueh the relevant protons 
(H) exehange with deuterium (D) in a solvent. 


Pr - H ± n 


solvent 


Pr - D ± H 


solvent 


( 1 ) 


Here Pr-H denotes a protein with a proton in the relevant hydrogen bond. The fraetionation ratio 
ean also be determined from UV speetroseopy ll^ . 


The fraetionation ratio is equilibrium eonstant of Eq. ([T]): 

^ ^ [Pr - D] [Hsolvent] 


( 2 ) 


[Pr - H] [Dsolvent] 

Translated into partition funetions, $ is essentially determined by the relative zero-point energy 


(ZPE) of a D relative to an H in the protein. As deseribed by Kreevoy and Eiang [Il3h . the ratio is 
given by 

fcfiT In $ = Zu-Pr — ^D-Pr ± -^D,solvent “ solvent (3) 

where T is the temperature and Zn-pr denotes the zero-point energy of a proton partieipating in 
the relevant hydrogen bond in the protein. Throughout this paper we set T = 300 K. 

Eraetionation is a purely quantum effeet. If the nuelear dynamies were elassieal, the frae¬ 
tionation ratio would be one. It would also be one if there were no ehanges in the vibrational 
frequeneies — more eorreetly, zero-point energies — of both H and D when they moved from the 
solvent to the protein. However, the vibrational potentials are different in the two environments. 
The donor-aeeeptor distanee, R, is typieally shorter in the protein, indieating a stronger H-bond 
and a softer X-H streteh potential (X is the H-bond donor). Consequently, the differenee between 
H and D zero-point energies gets larger llldll . and 4) gets smaller with deereasing R. However, 













3 


for very short bonds, typically when the donor and acceptor share the H or D atoms, the stretch 
frequencies begin to harden and $ starts to increase. $ then has a non-monotonic dependence on 

Riii. 


Mildvan and collaborators [[11 , 


1611 considered a particular parametrisation of the H-bond po¬ 


tential to connect the observed fractionation ratio with donor-acceptor bond lengths in a range of 
proteins. They generally find reasonable agreement between determinations of the length from the 
fractionation factor and that from the NMR chemical shift. In particular, the uncertainty is less 
than that deduced from X-rays. 

In this paper, we systematically investigate how the fractionation factor $ varies with the donor- 
acceptor distance R. Specifically, we consider the relative importance of different contributions to 
$. We find that the competing quantum effects associated with the X-H stretch and bend modes 
are the most significant. Non-degeneracy of the two bend modes and tunnel splitting of the stretch 
mode have small but noticeable effects. The only important effect of the secondary geometric 
isotope effect is that it enhances the contribution from the tunnel splitting, mostly for R ~ 2.4—2.6 
A. For most values of R, the value of we calculate differs from the empirical relation due to 
Mildvan et al. I116n that has been used to determine bond lengths in enzymes. 


II. METHOD 


We calculate the H/D zero-point energies, and hence $, with the electronic ground state poten¬ 
tial of a two-diabatic state model I117n . For X and Y as donor and acceptor, the two diabatic states 
of the model are X-H- ■ ■ Y and X- ■ ■ H-Y, which are modelled as Morse oscillators. The coupling 
between the diabats is a function of R, the X-Y distance, as well as the H-X-Y and H-Y-X angles; 
it decreases exponentially with increasing R and gets weaker with larger angular excursions of the 
H atom. Previo usly , we showed that this model can give a quantitative description of the correla¬ 
tions observed [Il8n for a diverse range of chemical compounds between R and X-H bond lengths, 
vibrational frequencies, and isotope effects [Il4ll . 

We now briefly discuss the domain of applicability of this simple model to hydrogen bonds 
in proteins, which are certainly complex and diverse chemical systems. First, our focus is on 
a small (but potentially important) sub-class of hydrogen bonds: short strong bonds. Second, 
we consider the simplest possible model that might capture the essential features of these bonds, 
independent of the finer structural details of a specific protein. The goal is to obtain physical insight 
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into the different quantum effeets that eontribute to the fraetionation factor, as well as their (non¬ 
monotonic) trends with donor-acceptor distance. H-bonds in proteins vary from weak to strong, 


and can further be modified by coupling to other neighbouring H-bonds ll^ . Also important 
are the proximity to and accessibility to the solvent, and anisotropic electric fields arising from 
neighbouring charged amino acid residues. An example of the latter occurs in the Photoactive 
Yellow Protein (PYP) where the existence of a possible low barrier H-bond may be dependent on 
deprotonation of the neighbouring Arg52 residue |l7i|201]. A key feature of the two-diabatic state 
model used here is that it takes the donor-acceptor bond distance and the pK^ difference as the key 
bond descriptors. These are input from available experimental information. These two parameters 
will certainly be modified by chemical substituents Kill , solvent, and perturbations from the local 
electric field as indicated above. Description of multiple H-bonds requires generalisation of the 
model considered here to include more than two diabatic states [l22n . 

The parametrization used in References [iJ and [l7| was for 0-H- • ■ O symmetric hydrogen 
bonds, i.e. the donor and acceptor have the same proton affinity (pK^). In the present work, 
we retain this parametrization. This is an approximation for comparisons with H-bonds in pro¬ 
teins, which are generally asymmetric (donor and acceptor with different pK^). Many H-bonds 
in proteins are actually N-H- ■ ■ O bonds. However, as the H-bonds become stronger (R < 2.5 A) 
the equal proton affinity approximation becomes more reliable. At such distances, the donor and 
acceptor effectively share the H atom. In the diabatic state model, the off-diagonal coupling ele¬ 
ment becomes large enough to strongly suppress or eliminate the barrier for the H atom transfer. 
Kreevoy and Liang lll3n . Bao et al. KSll . and Oltrogge and Boxer OlIl considered how asymmetry 
in the one-dimensional proton transfer potential modifies the fractionation factor. Non-degeneracy 
smaller than 800 cm“^ (or equivalently, a pK^ difference of about 2) has only a small effect on the 
proton transfer potential and the fractionation factor when R < 2.5 A. 

The total vibrational zero-point energy for Pr-H/D is 


Z{R) = Z||(i?) + Zx,o(i?) + (4) 

The three terms are the zero-point energy associated with X-H vibrations parallel to the hydro¬ 
gen bond (stretch), out-of-plane bend (o), and in-plane bend (i) of X-H- • ■ Y, respectively. (The 
plane typically refers to that of X-H.) The simple summation in the above equation points to our 
assumption that these modes are uncoupled. 

The 0-H/D stretch zero-point energy is calculated numerically by using the sine-function Dis- 














5 


Crete Variable Representation (DVR) [|24|] to solve the one-dimensional Sehrodinger equation as 


a funetion of R. This gives an essentially exaet treatment for the signifieant anharmonie and tun¬ 
nelling effeets that oeeur for low-barrier bonds [Il4l] . 

We treat zero-point energies of bend modes as half their elassieal harmonie frequeneies as a 
funetion of R. To break the degeneraey of the two modes, we use the result from the model itself 
that hardening of the two bend motions is similar: 

n^,o/^{Rf=u:l^o/i + V{R) (5) 

where oj \ „/,■ is the frequeney in the absenee of an H-bond and the funetion f{R) is given in Eqn. 
(6) of Ref. ilVi At least in the R range of interest, f{R) is a positive funetion that monotonieally 
deereases with inereasing R (eompare Figure 1). In general and so Dj_ > Here 

we take = 650 em“^ and = 650 or 1600 om“^. The uj parameters for the deuterium 
isotope are taken to be \/2 smaller than for the H isotope. 

Eq. ([3]) for $ assumes that only the ground state energies of the speeies is relevant (at the 
temperature of interest). For a symmetrie proton transfer potential, one expeets a first exeited state 
due to tunnel splitting that would be elose in energy to the ground state. In our model, this appears 
along the X-H streteh eoordinate. Henee the first exeited state for the H/D motion makes a further 
multiplieative eontribution of the form 

1 -f exp(-(5E'D/A:BT) 


$tun = 


( 6 ) 


1 -f exp(-5Ei^/fcBT) 
to the fraetionation faetor, where 6Eh/d = Eq- — Eq+ is the tunnel splitting. 

Another eontribution to $ eomes from the seeondary geometrie isotope effeet (SGIE) where the 
X-Y distanee ehanges upon deuteration. This is a subtle effeet with a non-monotonie dependenee 


on R [114 




2611 . As shown in Ref. [iJ and referenees therein, it arises beeause the rates of ehange 


with R of the zero-point energy for H and D are different (eompare Figure 1). The net effeet is 
that true minima of the total system energy with respeet to R for both H and D {Rh and Rd) 
are shifted relative to the elassieal minimum (Ro)- The differenee between the minima is largely 
under 0.04 A with Rd > Rh, but the resulting effeets on frequeneies are substantial. There are 
two eonsequenees of SGIE for the fraetionation ratio. First, the zero point energy for the H and 
D should be ealeulated at their respeetive minima. Seeond, an elastie energy assoeiated with the 
stretehing of the donor-aeeeptor distanee, of the form ^K{Rh/d — Rn)\ must be ineluded. Here, 
K is the elastie eonstant that is parametrised empirieally in Ref. [l4 We have ineluded both 
eonsequenees of the SGIE in our ealeulation. 
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The above details deseribe the ealeulation of the [Pr-D]/[Pr-H] part of in Eq. dH), as a funetion 
of R. The eorresponding ratio for the solvent is taken to be the ealeulated model value at i? = 2.8 
A, approximately the relevant length in water. Later we diseuss how our results are not partieularly 
sensitive to this exaet ehoiee of a referenee distanee. 

III. RESULTS AND DISCUSSION 



Figure 1. Competing quantum effects. The zero-point energies of the out-of-plane bend mode and of the 
two lowest stretch quantum states (due to tunnel splitting, Eq±) are shown for both H and D isotopes as 
a function of the donor-acceptor distance R. The black curves [H-stretch {Eq± )] correspond to the stretch 
mode of the H isotope, while green curves [D-stretch (Egi)] are those of the D isotope (solid: -|-, dot- 
dashed: —). Blue (H-bend) and magenta (D-bend) curves are out of plane bend zero-point energies for the 
H and D isotopes, respectively. With increasing R, the stretch zero-point energies increase and those of the 
bend decrease. Note how for R < 2.55 A (2.45 A), the tunnel splitting of the stretch mode for the H (D) 
isotope becomes observable. 


A. Role of competing quantum effects 

Figured! shows the computed trends of the stretch and bend zero-point energies from the model 
for both H and D isotopes. The solid black (H) and green (D) curves show the 0-H/D stretch 
zero-point energies. The limiting energies at large R are those of 0-H and 0-D bonds, about 1800 
and 1300 cm“^, respectively. Of relevance to the fractionation factor $ is the difference between 
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Figure 2. Effect of the bend modes and tunnel splitting on the fractionation ratio. The black curve (stretch) 
includes solely the effects of the X-H stretch vibrational mode. The blue curves include the effect of the 
in-plane (bj) and out-of-plane (bo) bending vibrational modes. The upper dashed curve (stretch, bo = bj) is 
for degenerate bend modes, while the lower solid curve (stretch, bo / bj, tunnel) includes the contribution 
from the first excited X-H stretch state (tunnel splitting) for non-degenerate bend modes. 


these eurves, whieh decreases (for the most part) with decreasing R. From Eq. ([3]), this leads to a 
drop in $. The black curve in Figure [2l shows the fractionation ratio with only the 0-H/D stretch 
zero-points included. 


A countervailing influence on $(i?) comes from the bends, which harden in frequency with 
decreasing R. Figure [U shows the zero-points only for the out-of-plane bends for the H and D 
cases. Their limiting values at large R are Iu!h,±,o (=325 cm”^) and -^ujh,±,o (=230 cm“^), 
respectively, where ujh,±,o = 650 cm~^. The corresponding trends for the in-plane bend (not 
plotted) are obtained from Eq. ([5]), setting = 1600 cm"k The consequences of the opposite 
trends for the hydrogen and deuterium bend and stretch frequencies — more compactly, competing 


quantum effects — has been the subject of much recent study lll4 . 


26 ^ 30 ]. 


Presently, for $, it is Zh — Zd that matters, which evidently also showcases the competi¬ 
tion between the X-H bends and the X-H stretch. The solid blue curve in Figure [2] shows that 
the hardening of the bend modes with decreasing R significantly increases the fractionation ratio 
compared to the contribution from just the stretch mode. This is one of the main results of this 


paper. Kreevoy and Fiang [1131] previously pointed out that bending modes could alter their results 
for the correlation of and R. They gave the rough estimate that $ could be increased by a factor 
of about 1.7. Edison, Weinhold, and Markley also mentioned the effect of the bend modes 113ill . 
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Figure 3. Corrections due to the secondary geometric isotope effect (SGIE) are shown as the dashed curves: 
the blue dashed line (stretch, bo / hi, SGIE) is without including tunneling and the red dashed line (stretch, 
bo / bj, tunnel, SGIE) is with tunneling. The solid curves of the same colour are from Eigure|2j which do 
not include the SGIE. The green curve is the empirical function, defined by Eqn. Q, and given in Ref. Il6. 
Note that the slight undulation in the dashed curves below 2.4 A can be traced to the very rapid change of 
the SGIE {Rd — Rh) in that region; see the lower panel of Eigure 7 in Ref. H. 


finding values of $ > 1 for weak bonds. 


B. Role of non-degeneracy of the bend modes 

For the main results of this paper, the bend modes are made non-degenerate using Eq. (|5]) since 
the out-of-plane and in-plane modes are different in frequency in reality. We briefly investigate 
how different the fractionation factor would be if the degeneracy were retained. 

The dotted blue curve in Figure [2] shows the plot of $ with both bends having frequencies 
corresponding to the out-of-plane mode. Evidently, the change is relatively small, but not neg¬ 
ligible. As per the model, the smaller the bend frequency, the faster it hardens; Vl±_{R)/u:±_ = 
Hence, changing one of the bends frequencies to 1600 cm reduces Zj^ 
relative to the degenerate case. Bend non-degeneracy can reduce the fractionation by about 10- 
20% compared to the degenerate case. 
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C. Role of the tunnel splitting 


For long H-bonds, the proton transfer 


potential has a high barrier and the tunnel splitting of 


the vibrational ground state is negligible lll4l] . However, for i? < 2.55 A, the splitting beeomes 


significant, as can be seen in Figure [B A multiplicative correction ^tun (Eq. ® ) introduces the 
effect of the thermal population of the first low-lying excited state. This factor is always larger than 
one because the H tunneling splitting is larger; SEu > 5 Ed- When the tunnel splitting is much 
larger than the thermal energy (i.e. for R < 2.4 A) the correction factor is extremely close 
to unity. When the tunnel splitting is much less than /c^T the correction factor is approximately 
(1 + {6En - dE^)/2kBT). 

Figured shows that the tunnel splitting has a small but non-negligible effect in the range, R ~ 
2.4 — 2.6 A. We should also clarify the nomenclature here. For sufficiently small R, the barrier 
is no longer present, and so there is strictly no “tunnel splitting.” We just have two well-separated 
vibrational energies instead. Note that in a solvent there will be local dynamical fluctuations of 
the local electric field that couple to the electric dipole moment associated with the X-H stretch 
and for large enough fluctuations the tunnel splitting will not appear because it will be destroyed 
by quantum decoherence [1320. Also, when the proton affinity of the donor and acceptor differ by 
more than about 500 cm“^ (1.5 kcal/mol or 1.3 pK^ units) this effect may be absent. 


D. Role of the secondary geometric isotope effect 


The SGIE has a significant effect on the stretch mode vibrational frequencies for i? ~ 2.4 — 
2.5 A, where the proton transfer barrier has effectively disappeared. Its inclusion yielded better 
agreement of the H/D stretch frequency ratio; compare Figure 8 in Ref. [l4|, where this ratio is 1 
(1.15) with (without) SGIE, which is a sizeable change for strong short H-bonds. 

EigurelH however, points to a more modest influence of the SGIE on the fractionation ratio. 
The only significant effect is how it modifies the contribution from the tunnel splitting. Without 
the SGIE, the correction factor is approximately (1 -f {6Eh{R) — 5Ed{R))/ 2,kBT). With the 
SGIE, the correction factor is approximately (1 -f {5Eh{Rh) — ^Ej:,{Rd)) /‘^^bT). This is larger 
because Rd > Rh, and an increase in donor-acceptor distance of as little as 0.02 A for D relative 
to H can increase the energy barrier and thereby noticeably decrease 5Ej:){Rd) llldll . 


If tunneling contributions are suppressed by, e.g., a sizeable difference in pKaS between the 
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donor and acceptor, the above analysis suggests that the SGIE would have only a small influence 
on the fraetionation factor. The key effect that appears to govern the magnitude range of $(i?) 
aeeording to our model is the competing quantum effeet between the X-H streteh and X-H bends. 


E. Sensitivity to choice of reference distance 

Calculation of <h, in Eqn. Q requires knowledge of the ratio of H/D eoneentration in the 
solvent. In the reported calculations we took this ratio to be given by the value ealeulated within 
our model at i? = 2.8 A, approximately the relevant length in water. Our results are not partieularly 
sensitive to this exact choice of this referenee distanee. Eor R > 2.7 A the differenee between the 
H and D zero-point energies is small. This can be seen from EigureO Even for the ease of purely 
streteh modes taking the reference distanee to be i? = 2.7 A would only inerease $(7?) by about 
10 per cent eompared to the values shown in our curves. 


F. Comparison with experiment 


To put our results in eontext we now briefly review previous measurements of $ that have been 
used to deduee a value for i? in a specific molecule. 


Based on ealeulations from an empirieal one-dimensional quartic potential [l23h . Mildvan et 


al. lllbn eonsidered a relation between the fraetionation faetor and donor-acceptor distance, 


R = (2.222 + 1.192$ - 1.335$^ + 0.608$^) A. 


(7) 


It was used together with measurements of fraetionation ratios for 18 H-bonds in several different 
enzymes to deduee the length R. The values they obtained for R were mostly in agreement with 
values of R dedueed from NMR chemical shifts, and from X-ray crystallography. Values of $ 
ranged from 0.32 to 0.97 and the eorresponding values of R were in the range 2.49 to 2.68 A. 
However, Eigure|3] shows significant differences between equation (|7]) and our results. 

Klug et al. [l33ll studied erystals of the dihydrated sodium salt of hydrogen bis(4-nitrophenoxide) 
and found a fraetionation ratio of 0.63 ± 0.04. Using the Kreevoy and Eiang I113n parametrisation, 
they noted that this value was ineonsistent with the bond length observed via X-rays, R = 2.452 A, 
and with the value of $ = 0.31 ±0.03 deduced from UV spectroseopy for bis(4-nitrophenoxide) in 
acetonitrile solution. Consequently, they suggested that “the solid was not in isotopic equilibrium 
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with the solvent from which it was precipitated.” However, their results are consistent with our 
parameterisation of $ versus R, if tunnel splitting is not included. 


Loh and Markley 11341] found fractionation factors in the range 0.28-1.47 for the different H- 


bonds in the protein staphylococcal nuclease. Clearly, we cannot explain their $ values larger 
than unity. However, it should be pointed out that many of these bonds are weak with donor- 
acceptor distances in the range, i?~2.8 — 3.3 A and that no clear correlation was observed 
between the values of and R. Loh and Markley did note a difference between the fractionation 
factors of backbone amide bonds that are solvent-exposed (average value 0.98) and those that are 
not (average value 0.79). This observation is consistent with what one might anticipate: solvent 
accessibilility could affect the the effective pK^, of the donor and/or acceptor. A consequent change 
(increase) in distance between them may also ensue. These would lead to a weakening of the H- 
bond, and therefore move values closer to 1. In the model we employ, this solvent effect would 
enter parametrically as a difference in effective pK^ and donor-acceptor distance. However, the 
works of Khare et al. OSH and LiWang and Bax 11361] point to minimal effects. The former, 
which reported on the immunoglobulin G binding domains of protein G, found little difference 
between average fractionation factors (average of 1.05 for a-helical, 1.13 for /3-sheet, and 1.08 
for solvent-exposed residues). LiWang and Bax gave similar findings for ubiquitin. A recent study 


oft 

£l. 


le core of protein Kinase A also found no correlation between $ values and secondary structure 


Thakur et al. 11381] have recently developed a new method for the rapid determination of H/D 


exchange from two-dimensional NMR spectra. Section S5 of their Supplementary material shows 
fractionation values for three different proteins. For 80 different amino acid residues in Tim23, 
the values ranged from 0.81 to 1.73. For 58 different amino acid residues in Ubiquitin, the values 
ranged from 0.34 to 1.67. For 54 different amino acid residues in Dph4, the values ranged from 
0.45 to 2.04. 


Recently, an extensive study was made of mutants of the Green F 


short H-bond between the chromophore and the amino acid Asp 148 112 ih . The donor-acceptor 


uorescent Protein with a 


bond length estimated from X-ray structures was 2.4 ± 0.2 A. The pK^ of the chromophore 
was systematically varied by 3.5 units through halogen substitutions. This range covers the pK^ 
matching (degenerate diabatic states) required for strongest bonds 11171] . The experimental results 
were compared to calculations based on a one-dimensional proton transfer potential based on same 
diabatic state model used here. The measured fractionation factors (deduced from analysis of UV 
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absorption spectra) were in the range 0.54 - 0.9, taking a minimum value for pK^ matching. This 
observation and a value of = 0.54 for R = 2.4 ± 0.2 A are consistent with our analysis when 
the bend modes and tunnel splitting are taken into account. 

Edison, Weinhold, and Markley performed ab initio calculations for a wide range of peptide 
clusters [13ill . They observed a correlation between the fractionation ratio and the donor-acceptor 
distance. For R > 2.55 A, the fractionation was larger than one, and for R ~ 2.45 A, $ ~ 0.6. 


IV. CONCLUSIONS 


We have shown that the H/D fractionation factor $ is quite sensitive to the donor-acceptor 
distance R in hydrogen bonds, and so, in principle, can be used as a “ruler” for determining 
bond lengths. However, caution is in order because there are a number of subtle effects that 
modify the exact form of the relationship between <h and R. These include competing quantum 
effects between stretch and bend modes, non-degeneracy of the bend modes, tunnel splitting, the 
secondary geometric isotope effect, and differences between the proton affinity of the donor and 
acceptor. 

Our results raise questions about whether values of $ as small as 0.3 are really possible for 
short bonds, contrary to some measurements and previous theoretical claims. Equally, our results 
cannot explain $ values that are much larger than 1 for long bonds. The discrepancy for short 
bonds may be due to our assumption that the stretch and bend modes are independent. Although 
our model quantitatively describes many experimental results for bond lengths, vibrational fre¬ 
quencies, and isotope effects, for R ~ 2.45 A, it does give stretch mode frequencies that are 
higher than observed. (See Figure 6 in Ref. [l4| and the associated discussion.) This would also 
lead to a larger fractionation factor than observed. Addressing this issue will require a system¬ 
atic investigation of solutions to the vibrational Schrodinger equation for a higher-dimensional 
(probably four dimensional) potential energy surface. We leave that for a future study. 
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